
global root_dir = "`1'"

include "$root_dir/code/config/config.do"


cap noi log using ${log_dir}/figure_4_macrosimulation.log, replace name(fig)

capture noi {

    global d = "`1'"
    do ${d}/code/macrosim/simulation_config.do "`1'/" 

    * Set figure specific configurations
    local mvers "medianp33"
    global N = 200

    * $N is the batch size for each simulation, but we have 4 batches (0-3), aggregated files are named according to total size
    global N = $N * 4  
    global lastyear = 2011

    * ----------------------------------------- *
    * Collapse runs
    * ----------------------------------------- *

    * append batch-combined models from simulation, go to wide

    use ${outdr}/data/runs_${ln_vers}pauto95_${chosen_spec}${FE}_recomputerecompute_0_${iwvers}${osmtrim}${spillrcvers}_n${N}.dta, clear
    ren pdepvar pdepvar_m0

    mmerge year run division using ${outdr}/data/runs_${ln_vers}pauto95_${chosen_spec}${FE}_recomputerecompute_1_${iwvers}${osmtrim}${spillrcvers}_n${N}.dta, unmatched(master)
    ren pdepvar pdepvar_m10

    mmerge year run division using ${outdk}/data/runs_${ln_vers}pauto95_${chosen_spec}${FE}_keepkeep_1_n${N}, unmatched(master)
    ren pdepvar pdepvar_m10_keepall

    mmerge year run division using ${outdk}/data/runs_${ln_vers}pauto95_${chosen_spec}${FE}_recomputekeep_1_n${N}, unmatched(master)
    ren pdepvar pdepvar_m10_keepspill

    mmerge year run division using ${outdk}/data/runs_${ln_vers}pauto95_${chosen_spec}${FE}_keepkeep_0_n${N}, unmatched(master) 
    ren pdepvar pdepvar_m0_keepall
    drop _m

    drop if year > ${lastyear}
    decode division, gen(sdiv)
    drop division
    sort year run sdiv

    reshape wide depvar pdepvar_*, i(year run) j(sdiv) string


    * calculate shares for mean, plus and minus 1 SD
    ren (depvarauto95 depvarpauto95) (auto95_data pauto95_data)
    gen share_data = auto95 / (auto95 + pauto95)
    foreach model in m0 m10 m10_keepall m10_keepspill m0_keepall {
        gen share_`model' = pdepvar_`model'auto95 / (pdepvar_`model'pauto95 + pdepvar_`model'auto95)
        bys year: egen share_`model'_sd = sd(share_`model')
        bys year: egen share_`model'_m = mean(share_`model')
        bys year: egen share_`model'_md = median(share_`model')
        gen share_`model'_u = share_`model'_m + share_`model'_sd
        gen share_`model'_l = share_`model'_m - share_`model'_sd
        gen share_`model'_mp5 = .
        gen share_`model'_mp10 = .
        gen share_`model'_mp25 = .
        gen share_`model'_mp50 = .
        gen share_`model'_mp90 = .
        gen share_`model'_mp95 = .
        gen share_`model'_mp75 = .
        gen share_`model'_mp33 = .
        gen share_`model'_mp66 = .
        forval year = 1997/2011 {
            su share_`model' if year == `year', d
            replace share_`model'_mp25 = `r(p25)' if year == `year'
            replace share_`model'_mp75 = `r(p75)' if year == `year'
            replace share_`model'_mp50 = `r(p50)' if year == `year'
            replace share_`model'_mp5 = `r(p5)' if year == `year'
            replace share_`model'_mp95 = `r(p95)' if year == `year'
            replace share_`model'_mp10 = `r(p10)' if year == `year'
            replace share_`model'_mp90 = `r(p90)' if year == `year'
            local p33
            local p66
            local third = 100*(1/3)
            local twothird = 100*(2/3)
            centile share_`model' if year == `year', centile(`third')
            local p33 = r(c_1)
            centile share_`model' if year == `year', centile(`twothird')
            local p66 = r(c_1)
            replace share_`model'_mp33 = `p33' if year == `year'
            replace share_`model'_mp66 = `p66' if year == `year'
        } 
    }

    * make a dummy for the nearest observation to the 5th percentile
    sort year run 
    su run
    local nruns = r(max) + 1

    xtile xt_share_m0 = share_m0 if year == 2011, nq(`nruns')
    foreach p in 5 25 33 50 66 75 95 { 
        gen dist_to_p`p' = abs(xt_share_m0 - `p')
        egen min_dist_to_p`p' = min(dist_to_p`p')
        replace min_dist_to_p`p' = . if min_dist_to_p`p' != dist_to_p`p'
        bys min_dist_to_p`p': gen p`p'auto95_ = 1 if _n == 1 & min_dist_to_p`p' != .
        bys run: egen p`p'auto95 = max(p`p'auto95_)
        replace p`p'auto95 = 0 if p`p'auto95 == .
        drop dist_to_p`p' min_dist_to_p`p' p`p'auto95_
    }
    
    * trace the p5, p25, p50, p75, p95 lines
    foreach p in 5 25 33 50 66 75 95 {
        foreach model in m0 m10 m10_keepall m10_keepspill m0_keepall {
            gen share_`model'_p`p'_ = pdepvar_`model'auto95 / (pdepvar_`model'pauto95 + pdepvar_`model'auto95) if p`p'auto95 == 1
            bys year: egen share_`model'_p`p' = max(share_`model'_p`p'_)
            drop share_`model'_p`p'_
        }
    }

    collapse (mean) share_*_* *data*, by(year)
    foreach var of varlist share_* {
        local var : subinstr local var "share_" "", all
        gen pct_`var' = 100 * share_`var'
    }

    * ----------------------------------------- *
    * Produce and export figure
    * ----------------------------------------- *

    include ${code_dir}/config/figuretools.do

    qui colorpalette scico bamO, n(6) saturate(40) intensity(0.9)
    di "`r(p1)'"
    di "`r(p2)'"
    di "`r(p3)'"
    di "`r(p4)'"
    di "`r(p5)'"
    di "`r(p6)'"

    local v pct

    local vname "(%)"
    local range 6 20
    local textpos 8
    local mvers_dis "_`mvers'"
    
    local diagntext "text(`textpos' 2008 "Wagechange: $wc. Spec: $chosen_spec$FE" "Spill: ${iwvers}${osmtrim}${spillrcvers}" "N keep: ${N}" "Shaded areas are +/- 1 sd", size(vsmall) just(left))"
    
    local u mp33
    local l mp66
    local m mp50
    
    twoway (rarea `v'_m0_`u' `v'_m0_`l' year, color("150 193 74*.9%50") lw(none)) /// + baseline 
        || (rarea `v'_m10_keepall_`u' `v'_m10_keepall_`l' year, color("237 96 208*.9%50") lw(none)) /// direct 
        || (rarea `v'_m10_`u' `v'_m10_`l' year, color("87 117 0*.9%40") lw(none)) /// total
        || (rarea `v'_m10_keepspill_`u' `v'_m10_keepspill_`l' year, color("255 170 187*.9%50") lw(none)) /// direct + stock
        || (connected `v'_data `v'_m0_`m' `v'_m10_keepall_`m' `v'_m10_keepspill_`m' `v'_m10_`m' year, ///
                msymbol(d o d d o) msize(small small small small small small) ///
                lwidth(thin thin thin thin thin thin) lpattern(shortdash solid solid solid solid solid solid) ///
                color("`r(p1)'" "150 193 74*.9" "237 96 208*.9"  "255 170 187*.9" "87 117 0*.9") ///
                xtitle("") xlabel(1998(3)2010, nogrid) ///
                ytitle("Automation `vname'") ylabel(0(5)25, grid glpattern(shortdash) glwidth(thin) glcolor(gs15)) yscale(range(`range')) /// 
                legend(order(5 "Data" 6 "Baseline" 7 "+ direct" 9 "+ total" 8 "+ stock") ///
                    ring(2) pos(6) size(small) cols(3) region(lstyle(solid) lw(none) lcol(gs8))) ///
                plotregion(margin(0))) 
    graph export ${fig_dir}/main/Figure_4_macrosimulation_${ln_vers}pauto95_${chosen_spec}${FE}_${iwvers}${osmtrim}${spillrcvers}_${wc}_n${nr}`v'`mvers_dis'`drop_outliers_dis'.pdf, replace
    graph export ${fig_dir}/main/Figure_4_macrosimulation_${ln_vers}pauto95_${chosen_spec}${FE}_${iwvers}${osmtrim}${spillrcvers}_${wc}_n${nr}`v'`mvers_dis'`drop_outliers_dis'.eps, replace

    cap log using ${log_dir}/figure_4_macrosimulation_numbers.log, replace name(fig_numbers)

    *reporting some numbers
    keep if year > 1996 & year < 2012
    gen total_baseline_gap = `v'_m10_`m' - `v'_m0_`m'
    summarize total_baseline_gap
    di "mean total_baseline_gap: `r(mean)'"
    gen direct_data_gap = `v'_m10_keepall_`m' - `v'_data
    summarize direct_data_gap
    di "mean data_direct_gap: `r(mean)'"
    summarize `v'_data if year == 1997
    di "mean share_data 1997: `r(mean)'"
    summarize `v'_data if year == 2011
    di "mean share_data 2011: `r(mean)'"
    cap log close fig_numbers
}
if _rc == 0 {
    display "Execution finished successfully."
}
else {
    display "Execution finished with errors."
}

cap log close fig